Aquest pipeline s’ha dissenyat per tal de poder realitzar un anàlisi integratiu de dades multiòmiques, partint de dades de proteòmica i de transcriptòmica. Tot i així, podria ser adaptat per a incloure-hi dades de metabolòmica o RNA-seq. Abans de començar l’anàlisi, cal estar segur que les mostres que s’inclouen als dos anàlisi coincideixen i estan en la mateixa codificació.

L’anàlisi consta de 3 aproximacions diferents. En la primera aproximació, es comparen els gens i proteïnes diferencialment expressades (amb el cut-off escollit) per tal d’identificar-ne de comuns. En la segona aproximació, es realitza un anàlisi multivariant anomenat multiple co-inertia analysis (MCIA) (Meng, 2014) per tal d’integrar múltiples capes de dades -òmiques. Per últim, també es realitza de manera complementària un anàlisi d’enriquiment combinat, realitzant un multiGSEA (Canzler, 2020).

Per a l’anàlisi multiòmic, en primer lloc es necessita tenir els gens i les proteïnes en la mateixa codificació. Per aquest motiu, transformarem la codificació Uniprot de les proteïnes a Gene Symbol, amb el paquet AnnotationDbi (Pagès 2019).

prots_all <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/TT_MCAO_all_prot.csv", header=TRUE, sep=",")
prots_all$X <- substr(prots_all$X, 1, 6)  
rownames(prots_all) <- prots_all$X
library(AnnotationDbi)
library(org.Mm.eg.db)
AnnotationDbi::keytypes(org.Mm.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GO"           "GOALL"        "IPI"          "MGI"          "ONTOLOGY"    
## [16] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         "PROSITE"     
## [21] "REFSEQ"       "SYMBOL"       "UNIGENE"      "UNIPROT"
symbols_prot <- AnnotationDbi::mapIds(org.Mm.eg.db,
                                            keys = rownames(prots_all), keytype = "UNIPROT",
                                            column = c("SYMBOL"))
na.symbols_a <- which(is.na(symbols_prot))
prots_symbol_all <- cbind(prots_all, symbols_prot)

write.csv(prots_symbol_all, "E:/Dades crues Multi-Bio-Targets/Integromics/TT_proteomics_symbol_all.csv")

S’hauran d’anotar manualment aquelles proteïnes que no s’han pogut recodificar amb el paquet. Per això, s’utilitza la web d’ Uniprot.

Atenció! En aquest anàlisi estem treballant amb mostres de ratolí. Si estem treballant amb mostres d’una altra espècie, com per exemple humà, haurem d’utilitzar les anotacions humanes (org.Hs.eg.db). Això s’ha de tenir en compte també a l’hora de fer l’anotació manual amb Uniprot.

Comparació de molècules comuns

Ara passem a comparar les dues llistes de molècules diferencialment expressades. Per això, s’ha també anotat la llista de proteïnes significatives, i posteriorment aquelles que no s’han pogut anotar s’ha fet manualment com en el cas anterior.

prots <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/TT final proteomics_int.csv", header=TRUE, sep=";")
rownames(prots) <- prots$X
AnnotationDbi::keytypes(org.Mm.eg.db)
##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
##  [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
## [11] "GO"           "GOALL"        "IPI"          "MGI"          "ONTOLOGY"    
## [16] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         "PROSITE"     
## [21] "REFSEQ"       "SYMBOL"       "UNIGENE"      "UNIPROT"
symbols_prots <- AnnotationDbi::mapIds(org.Mm.eg.db,
                                            keys = rownames(prots), keytype = "UNIPROT",
                                            column = c("SYMBOL"))
na.symbols <- which(is.na(symbols_prots))
prots_symbol <- cbind(prots, symbols_prots)

write.csv(prots_symbol, "E:/Dades crues Multi-Bio-Targets/Integromics/TT_proteomics_symbol.csv")

En aquest pas és on es comparem les dues llistes, un cop tenen la mateixa codificació.

prots_symbol_man <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/TT_proteomics_symbol_manually_sign.csv", header=TRUE, sep=";")
genes_symbol <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/TT final male_gens_int.csv", header=TRUE, sep=";")

protsymb <- prots_symbol_man$symbolsFromUniprot
genesymb <- genes_symbol$SYMBOL
common_molecules <- intersect(protsymb, genesymb)
length(common_molecules)
## [1] 0

En el cas d’aquests dos datasets, no hi ha molècules significatives comuns. Això podria ser degut a que no hi ha gaire anàlits estudiats coincidents. Ho comprovem amb un diagrama de Venn.

library(VennDiagram)
base_01 <- read.delim("E:/Dades crues Multi-Bio-Targets/venn diagrams integromics/comparison.txt", header = TRUE, sep = "\t")
A <- base_01$SYMBOL
B <- base_01$symbols_prot
B <- B[-c(2486:5071)]

a1=length(A)
b1=length(B)

ab=intersect(A, B)
ab1=length(ab)
draw.pairwise.venn(a1, b1, ab1, fill = c("skyblue", "pink1"), lty = "blank", title="raw p-value<0.05")
Diagrama de Venn de les molècules analitzades en cadascun dels casos, després dels respectius filtratges

Diagrama de Venn de les molècules analitzades en cadascun dels casos, després dels respectius filtratges 

## (polygon[GRID.polygon.1], polygon[GRID.polygon.2], polygon[GRID.polygon.3], polygon[GRID.polygon.4], text[GRID.text.5], text[GRID.text.6], text[GRID.text.7], lines[GRID.lines.8], text[GRID.text.9], text[GRID.text.10])

Com veiem al diagrama, només hi ha 179 molècules coincidents en els dos anàlisis, la qual cosa podria explicar perquè no en trobem de significatives comuns.

MCIA

Per a l’anàlisi multivariant, s’ha escollit realitzar un anàlisi de co-inèrcia múltiple o MCIA (Meng, 2014), que consisteix en un mètode exploratori d’anàlisi de dades per tal d’identificar co-relacions entre diferents datasets multidimensionals, integrant diferents capes d’informació. Per a realitzar-lo, necessitem pujar les matrius de les dades normalitzades dels dos datasets, assegurant-nos que les mostres estan ordenades de la mateixa manera, i que les variables tenen la mateixa codificació.

#BiocManager::install("omicade4")
library(omicade4)

GENES_mat <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/GENS_matrix.csv", header=TRUE, sep=",")
rownames(GENES_mat) <- GENES_mat$SYMBOL
GENES_mat <- GENES_mat[,-1]
PROTS_mat <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/proteomics.csv", header=TRUE, sep=",")
rownames(PROTS_mat) <- PROTS_mat$symbols_prot
PROTS_mat <- PROTS_mat[,-1]

data.list <- list(GENES_mat, PROTS_mat)

sapply(data.list, dim)
##      [,1] [,2]
## [1,] 5071 2485
## [2,]   14   14
mcia_male <- mcia(data.list, cia.nf = 2, cia.scan = FALSE, nsc = T, svd = TRUE)

Després de realitzar el MCIA, obtenim els gràfics resultants.

plot(mcia_male, axes = 1:2,
     sample.lab = TRUE, sample.legend = TRUE, sample.color = 2,
     phenovec = NULL, df.color = 2,
     df.pch = NA, gene.nlab = 5)
Anàlisi de MCIA

Anàlisi de MCIA

En el primer panell (superior esquerra), es presenta una projecció del MCIA on es representen els dos datasets (df1 és el dataset de transcriptòmica i df2 és el dataset de proteòmica). La llargada de les línies representa la divergència entre les mostres. Quant més curta és la línia, més concordança entre elles. En el cas d’exemple, la distància és similar en la majoria de les mostres, i no s’observen unes mostres molt més divergents que d’altres. Es pot observar una agrupació entre els individus 1,2,3 i 5 (esquerra) i una altra entre els individus 6, 7 i 8 per part del primer component

En el panell superior dreta, es representa quines variables dels dos datasets contribueixen més a la variància. Aquelles amb el major pes són les que es troben més allunyades del centre del gràfic, i es pot observar que dels que més aporten, que són els que estan etiquetats.

En el panell inferior esquerra, es representa el pes de cada un dels components sobre la variància total. L’eix de l’esquerra són els eigen values, i l’eix de la dreta representa el %.

Per últim, en el panell inferior dreta estan representats els pseudoeigen values, per tal d’indicar la similitud entre datasets i veure quin dels dos aporta més a la variància. Es pot observar que els dos datasets no són gaire similars.

A continuació, s’han seleccionat les 100 variables que més aporten de cadascun dels eixos del primer component (positius i negatius). Se n’han seleccionat 100 per cada eix i per cada data set (400 en total). Amb això, s’ha realitzat un over-representation analysis (ORA). Seleccionant aquestes variables, podem ser capaços de tenir una millor cobertura de l’anàlisi d’enriquiment, que no pas si escollíssim individualment de cada dataset.

select_mol100 <- topVar(mcia_male, axis=1, end="both", topN=100)
write.csv(select_mol100, "E:/Dades crues Multi-Bio-Targets/Integromics/selected_molecules_mcia100.csv")

S’utilitza el paquet d’R del Gprofiler (Kolberg, 2020) per fer l’ORA. En primer lloc hem de pujar el background de l’anàlisi, és a dir el total d’anàlits estudiats.

background <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/background_nodupl.csv", header=FALSE, sep=",")
selected_mol_column_100 <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/selected_molecules_mcia100_column.csv", header=FALSE, sep=",")

Posteriorment es duu a terme l’anàlisi, seleccionant les bases de dades de vies que volguem triar (en aquest cas la secció Biological Processes del Gene Ontology (GO) (Ashburner, 2000) i KEGG (Kanehisa, 2000)). Definim també l’organisme amb el qual treballem, el mètode de correcció que volem aplicar i el threshold.

library(gprofiler2)
background_vec <- c(background$V1)
selected_mol_column_100_vec <- selected_mol_column_100$V1
gprof_100 <- gost(selected_mol_column_100_vec, organism="mmusculus", correction_method="fdr", domain_scope="custom", custom_bg=background_vec, source=c("GO:BP", "KEGG"), user_threshold = 0.5)

Per tal de visualitzar gràficament els resultats, podem construir un Manhattan plot. En el cas de les nostres dades, només trobem 1 via enriquida després de l’anàlisi. Es tracta d’un gràfic interactiu, on clicant a cadascun dels cercles es pot veure de quina via es tracta.

gostplot(gprof_100)

Manhattan plot de l’anàlisi ORA

multiGSEA

Per tal de complementar l’anàlisi d’enriquiment anterior i utilitzar una altra aproximació, el que es realitza és un anàlisi GSEA ( gene set enrichment analysis ) combinat (Canzler, 2020). Típicament, en el GSEA es s’utilitza un llistat de tots els anàlits estudiats ordenat per algun estadístic (logFC, estadístic T, etc.). L’elecció de l’estadístic està a mans del analista, tot i que Zyla et al. (Zyla, 2017) van publicar un estudi en el qual presentaven l’impacte que pot tenir la tria de l’estadístic en el resultat final.

En primer lloc, necessitarem la llibreria de l’espècie que estem utilitzant.

#BiocManager::install(version='devel')
#BiocManager::install("multiGSEA")
library(org.Mm.eg.db)
library(multiGSEA)
library(magrittr)

A continuació, es pugen les llistes rankejades per separat. Cada arxiu haurà de tenir 3 columnes, una amb el ID del gen/proteïna, una altra amb el logFC i una úlima amb el p-valor.

transcriptome <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/multigsea_transcriptomics.csv", header=TRUE, sep=",")
proteome <- read.csv("E:/Dades crues Multi-Bio-Targets/Integromics/multigsea_proteomics.csv", header=TRUE, sep=",")

Amb la funció rankFeatures, creem l’estadístic ls, el qual té la següent fórmula: ls=sign(log2(FC))∗log10(p−value).

omics_data <- initOmicsDataStructure( layer = c("transcriptome", 
                                                "proteome"))
omics_data$transcriptome <- rankFeatures(transcriptome$logFC, 
                                         transcriptome$P.Value)
names(omics_data$transcriptome) <- transcriptome$SYMBOL
omics_data$proteome <- rankFeatures(proteome$logFC, proteome$P.Value)
names(omics_data$proteome) <- proteome$SYMBOL

Es mostren els 5 primers elements de cadascun dels llistats.

omics_short <- lapply(names(omics_data), function(name){ 
  head(omics_data[[name]])
})
names(omics_short) <- names(omics_data)
omics_short
## $transcriptome
##     Npas4      Ccl3       Fos      Fosb      Ccn1      Atf3 
## 10.149967 14.761954 10.903090 11.109579 10.958607  8.962574 
## 
## $proteome
##      Txn2       Hba    Hbb-b1       Sfn  Map1lc3a    Ppp1r2 
## 1.5174278 1.5222740 1.7768823 0.5734391 1.0659581 0.5107317

Posteriorment, necessitem descarregar les bases de dades que utilitzarem. Es basa en el paquet graphite (Sales, 2012), el qual té disponible les següents bases de dades:

Per aquest anàlisi, seleccionarem KEGG (Kaneisha, 2000), PathBank (Wishart, 2020) i Reactome (Jassal, 2020).

databases <- c( "kegg", "reactome", "pathbank")
layers <- names(omics_data)
pathways <- getMultiOmicsFeatures(dbs=databases,
                                  layer = layers,
                                   returnTranscriptome = "SYMBOL",
                                   returnProteome = "SYMBOL",
                                  organism="mmusculus",
                                   useLocal = FALSE)

Un cop tenim les bases descarregades, realitzem el GSEA per cada capa de dades per separat.

enrichment_scores <- multiGSEA(pathways, omics_data)

En aquest últim apartat, es calcularan els p-valors agregats, amb la funció extractPvalues. El que fa aquesta funció és crear un dataframe on cada fila representa una via i les columnes representen els p-valors i els p-valors ajustats de cada dataset. Així, es calcula el p-valor agregat i posteriorment el p-valor ajustat.

Per defecte, el mètode per calcular els p-valors agregats és el mètode Z o de Stouffer (Stouffer, 1949), però hi ha altres mètodes com el de Fisher (Fisher, 1932) o el de Edgington (Edgington, 1972). Un resum d’aquests mètodes es pot trobar a la publicació original de multiGSEA (Canzler, 2020).

df <- extractPvalues(enrichmentScores = enrichment_scores,
                      pathwayNames = names(pathways[[1]]))

df$combined_pval <- combinePvalues(df)
df$combined_padj <- p.adjust(df$combined_pval, method = "BH")

df <- cbind(data.frame( pathway = names(pathways[[1]])), df)

Ara ja podem extreure la taula amb els p-valors individuals de cada dataset i els p-valors agregats. Es representen les primeres files de la taula de resultats.

df2 <- df[order(df$combined_pval),]
write.csv(df2, "E:/Dades crues Multi-Bio-Targets/Integromics/resultats multiGSEA.csv")
df3 <- df2[1:10,]
knitr::kable(df3, booktabs = TRUE, caption = "Primeres files de l'anàlisi combinat multiGSEA")
Primeres files de l’anàlisi combinat multiGSEA
pathway transcriptome_pval transcriptome_padj proteome_pval proteome_padj combined_pval combined_padj
98 (KEGG) MAPK signaling pathway 0.0000005 0.0057058 0.7885986 0.9327148 0.2324486 0.9829597
168 (KEGG) IL-17 signaling pathway 0.0012642 0.3051784 NA NA 0.3051784 0.9829597
179 (KEGG) Circadian entrainment 0.0000114 0.0395878 0.2038664 0.9327148 0.4272968 0.9829597
505 (REACTOME) Signaling Pathways 0.0000093 0.0395878 0.2054264 0.9327148 0.4272968 0.9829597
301 (KEGG) Rheumatoid arthritis 0.0028212 0.5631564 NA NA 0.5631564 0.9829597
161 (KEGG) Toll-like receptor signaling pathway 0.0000880 0.1141167 0.7616387 0.9327148 0.5816238 0.9829597
244 (KEGG) Amphetamine addiction 0.0000839 0.1141167 0.0072102 0.9327148 0.5816238 0.9829597
543 (REACTOME) Signaling by NTRKs 0.0000839 0.1141167 0.8590604 0.9327148 0.5816238 0.9829597
1101 (REACTOME) MAPK targets/ Nuclear events mediated by MAP kinases 0.0000566 0.1141167 0.6493776 0.9327148 0.5816238 0.9829597
104 (KEGG) cAMP signaling pathway 0.0001120 0.1162175 0.3941980 0.9327148 0.5846075 0.9829597

En el cas de les dades d’exemple, cap via té un p-valor ajustat combinat menor a 0.05.

Informació de la sessió

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.1252    
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] magrittr_2.0.1       multiGSEA_1.0.2      gprofiler2_0.2.0    
##  [4] omicade4_1.30.0      ade4_1.7-16          VennDiagram_1.6.20  
##  [7] futile.logger_1.4.3  org.Mm.eg.db_3.12.0  AnnotationDbi_1.52.0
## [10] IRanges_2.24.1       S4Vectors_0.28.1     Biobase_2.50.0      
## [13] BiocGenerics_0.36.1 
## 
## loaded via a namespace (and not attached):
##   [1] TH.data_1.0-10              fgsea_1.16.0               
##   [3] colorspace_2.0-1            ellipsis_0.3.2             
##   [5] XVector_0.30.0              GenomicRanges_1.42.0       
##   [7] bit64_4.0.5                 fansi_0.5.0                
##   [9] mvtnorm_1.1-1               mathjaxr_1.4-0             
##  [11] codetools_0.2-16            splines_4.0.2              
##  [13] mnormt_2.0.2                cachem_1.0.5               
##  [15] knitr_1.33                  TFisher_0.2.0              
##  [17] jsonlite_1.7.2              graph_1.68.0               
##  [19] graphite_1.36.0             compiler_4.0.2             
##  [21] httr_1.4.2                  backports_1.2.1            
##  [23] Matrix_1.2-18               fastmap_1.1.0              
##  [25] lazyeval_0.2.2              formatR_1.11               
##  [27] htmltools_0.5.1.1           prettyunits_1.1.1          
##  [29] tools_4.0.2                 gtable_0.3.0               
##  [31] glue_1.4.2                  GenomeInfoDbData_1.2.4     
##  [33] dplyr_1.0.6                 rappdirs_0.3.3             
##  [35] fastmatch_1.1-0             Rcpp_1.0.6                 
##  [37] vctrs_0.3.8                 multtest_2.46.0            
##  [39] crosstalk_1.1.1             made4_1.64.0               
##  [41] xfun_0.23                   stringr_1.4.0              
##  [43] rbibutils_2.1.1             lifecycle_1.0.0            
##  [45] gtools_3.8.2                zlibbioc_1.36.0            
##  [47] MASS_7.3-51.6               zoo_1.8-9                  
##  [49] scales_1.1.1                hms_1.1.0                  
##  [51] MatrixGenerics_1.2.1        SummarizedExperiment_1.20.0
##  [53] sandwich_3.0-1              lambda.r_1.2.4             
##  [55] RColorBrewer_1.1-2          yaml_2.2.1                 
##  [57] memoise_2.0.0               gridExtra_2.3              
##  [59] ggplot2_3.3.3               stringi_1.6.2              
##  [61] RSQLite_2.2.7               highr_0.9                  
##  [63] mutoss_0.1-12               plotrix_3.8-1              
##  [65] checkmate_2.0.0             caTools_1.18.2             
##  [67] BiocParallel_1.24.1         GenomeInfoDb_1.26.7        
##  [69] Rdpack_2.1.2                rlang_0.4.11               
##  [71] pkgconfig_2.0.3             matrixStats_0.59.0         
##  [73] bitops_1.0-7                evaluate_0.14              
##  [75] lattice_0.20-41             purrr_0.3.4                
##  [77] htmlwidgets_1.5.3           bit_4.0.4                  
##  [79] tidyselect_1.1.1            R6_2.5.0                   
##  [81] snow_0.4-3                  gplots_3.1.1               
##  [83] generics_0.1.0              multcomp_1.4-17            
##  [85] DelayedArray_0.16.3         DBI_1.1.1                  
##  [87] pillar_1.6.1                sn_2.0.0                   
##  [89] survival_3.1-12             scatterplot3d_0.3-41       
##  [91] RCurl_1.98-1.3              tibble_3.1.2               
##  [93] crayon_1.4.1                futile.options_1.0.1       
##  [95] KernSmooth_2.23-17          utf8_1.2.1                 
##  [97] tmvnsim_1.0-2               plotly_4.9.3               
##  [99] rmarkdown_2.8               progress_1.2.2             
## [101] data.table_1.14.0           blob_1.2.1                 
## [103] metap_1.4                   digest_0.6.27              
## [105] tidyr_1.1.3                 numDeriv_2016.8-1.1        
## [107] munsell_0.5.0               viridisLite_0.4.0

Bibliografia

Ashburner et al. 2000. Gene ontology: tool for the unification of biology. Nat Genet. 25(1):25-9.

Canzler, S. and Hackermüller, J. 2020. multiGSEA: a GSEA-based pathway enrichment analysis for multi-omics data. BMC bioinformatics, 21(1), 561.

Edgington, Eugene S. 1972. An Additive Method for Combining Probability Values from Independent Experiments. The Journal of Psychology 80 (2): 351–63.

Fisher, S R A. 1932. Statistical Methods for Research Workers - Revised and Enlarged. Edinburgh, London.

Jassal, B., Matthews, L., Viteri, G., Gong, C., et al. 2020. The reactome pathway knowledgebase. Nucleic Acids Res. 48(D1):D498-D503.

Kanehisa, M. and Goto, S. 2000. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28, 27-30.

Kolberg, L. and Raudvere, U. 2020. gprofiler2: Interface to the ‘g:Profiler’ Toolset. R package version 0.2.0. https://CRAN.R-project.org/package=gprofiler2

Meng, C., Kuster, B., Culhane, A. C., and Gholami, A. M. 2014. A multivariate approach to the integration of multi-omics datasets. BMC bioinformatics, 15, 162.

Pagès, H., Carlson, M., Falcon, S and Li, N. 2019. AnnotationDbi: Manipulation of Sqlite-Based Annotations in Bioconductor.

Sales, G., Calura, E., Cavalieri, D. and Romualdi, C. 2012. graphite - a Bioconductor package to convert pathway topology to gene network. BMC Bioinformatics, 13, 20.

Wishart, D. S., Li, C., Marcu, A., et al. 2020. PathBank: A Comprehensive Pathway Database for Model Organisms. Nucleic Acids Res. 48(D1):D470-D478.

Zyla, J., Weiner, J., Marczyk, M. and Polanska, J.. 2017. Ranking Metrics in Gene Set Enrichment Analysis: Do They Matter? BMC Bioinformatics 18 (256).